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The generalized master equation or the equivalent continuous time random walk equations can 
be used to compute the macroscopic first passage time distribution (FPTD) of a complex stochastic 
system from short-term microscopic simulation data. The computation of the mean first passage time 
and additional low-order FPTD moments can be simplified by directly relating the FPTD moment 
generating function to the moments of the local FPTD matrix. This relationship can be physically 
interpreted in terms of steady-state relaxation, an extension of steady-state flow. Moreover, it is 
amenable to a statistical error analysis that can be used to significantly increase computational 
efficiency. The efficiency improvement can be extended to the FPTD itself by modelling it using a 
Gamma distribution or rational function approximation to its Laplace transform. 
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I. INTRODUCTION 

The first passage time distribution (FPTD) Q] con- 
cisely describes the kinetics of macroscopic transitions of 
complex macromolecular systems; e.g., the transition of 
a disordered heteropolymer from random to specifically- 
absorbed conformations Q or the dynamics of protein 
folding 0. 0] . If a statistical ensemble of systems is pre- 
pared at time t = in an initial metastable macroscopic 
state (macrostate) i and Pf (t) is the probability that an 
ensemble member is in absorbing final macrostate / at 
time t, then the FPTD is 



<p(r) = dP f (t)/dt\ t=T 



(1) 



We assume that / is the only absorbing state and that 
the system is ergodic, so 



ip(r)dT — Pf(oo) = 1 



(2) 



The mean first passage time (MFPT) is the first moment 
(rip) = ((t)), where 



oc 

:./■) = / x(r) dr . 

10 

/>oo 

((x)) = / x(T)lfi(T)dT. 

Jo 



It determines the transition rate as ((r)) _1 [EE1- 

Except for the simplest models, (p(r) and its moments 
can not be analytically computed. Moreover, direct nu- 
merical computation, e.g. by molecular or stochastic dy- 
namics, is often unaffordable. For example, proteins typ- 
ically have 10 3 — 10 5 conformational degrees-of-freedom 



and the macroscopic timescales of interest can be > 10 12 
times larger than the microscopic timescale 0, so huge 
amounts of computational effort would be needed for di- 
rect simulation. 

Coarse-graining can overcome these problems. The 
essential idea is to subdivide the macroscopic transi- 
tion into a network of discrete intermediate mesoscopic 
transitions that are fast enough for feasible computation 
(e.g., using Monte Carlo @ or molecular dynamics 
methods) yet slow enough (relative to the microscopic 
timescale) for approximation by a first-order stochastic 
equation. 

The simplest approximation of this sort is a Markovian 
master equation [l| for P(t), the N- vector that specifies 
ensemble probability over the intermediate mesoscopic 
states and initial and final states (see [l(| and references 
therein for examples in the context of protein folding). 
However, this approach will only be accurate when each 
mesoscopic transition can be characterized as a simple 
Poisson process. This is not the case for many important 
problems because of fractal or quasi-diffusive dynamics or 
because the timescale that would be needed to achieve 
the Markovian limit is too long for direct computation 

mm. 

The non-Markovian classical generalized master equa- 
tion is a more robust app roximation that can be used in 
this situation Assuming injection at t = of 

ensemble members into initial state i, it is 



dP{t) 
dt 



5(t - (T 1 



p CO 

/ T(r) • P(t - r)dr , 
Jo 



(3) 
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where T(t) is the N x N matrix of transition functions 
which include memory effects and e s denotes the basis 
vector that has component s equal to 1 and all other 
components equal to 0. The first term, along with the 
boundary condition P{— oo) = 0, initializes the system 
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with no memory at t = 0. Conservation of probability 
and causality imply that 

l-r(r)=0, IV s (t)<0 (s'/s), 

(where 1 is the TV-vector with all components equal to 
1) and s and s' take values corresponding to i, /, and all 
the intermediate mesoscopic states. Therefore 



1 • dP(t)/dt = S(t - 0+) 
Since / is absorbing, 



r(r) 



= 



(4) 



(5) 



If T(t) can be analytically computed (e.g., by projection 
0)0)) then Eqs. and © can be used to compute 
<p{r). 

When r(r) can not be computed, an approach us- 
ing numerical simulations can be employed [9|. By ini- 
tializing multiple simulations in state s and determin- 
ing the distribution of waiting times for first-transitions 
to the other mesoscopic states, simulations can be used 
to determine K(t), the N x N local FPTD matrix 
(sometimes called the "first-jump waiting time" matrix). 
—K s i s (t) (s 7^ s') is the probability density that, af- 
ter arriving at state s, a system waits for an interval 
r before first leaving and that it goes to s'. For com- 
pact notation, we define the diagonal elements of K(t) 
as K ss {t) = - J2 S ^ S K s , s (t). Thus, like T(r), K(r) sat- 
isfies 

1-Jf(r) = 0, K ss ,(r)<0 {s' + s), K(t) ■ e f = , 

and, by its definition and the assumption that / is the 
only absorbing state, 



/>oo 

/ K ss {r)dT = l (sjtf) 
Jo 



(7) 



K (t) is the kernel of an alternative representation of 
stochastic dynamics with memory — the generalized con- 
tinuous time random walk (CTRW). Originally intro- 
duced to describe random walks on lattices [lg, [l9| , the 
CTRW was generalized |2jJ, EI to a form that can be 
extended to memory-dependent stochastic processes on 
a mesoscopic network of arbitrary connectivity. In our 
notation and assuming t = initialization in state i this 
is 



dP(t) 
dt 



Q(t) 



/ Ko{t) ■ Q(t — r)dr 
Jo 



poo 

Q(t) = <5(t-0+)e J + / Kp(r)-Q(jt-T)dT(jSb) 
Jo 

where Kd and Kp are the matrices comprised, respec- 
tively, of the diagonal or off-diagonal elements of K, and 
the boundary conditions are P{— oo) = Q(— oo) = 0. Eq. 
(|8b|) implies that Q s {t)dt is the probability that an en- 
semble member makes a transition to state s within the 



interval [t,t + dt), and Eq. (|8a|l states that dP(t)/dt is 
the difference between the incoming and outgoing prob- 
ability flows. 

Eqs. iJSJ and © provide equivalent descriptions of the 
temporal evolution of P(t) |13lll5l |. and comparing their 
Laplace transforms shows that 

f(u) = uK{u) -[I - Kd(u)]- 1 (9a) 
K(u) = T(u) ■ [ul + f D (u)}- 1 , (9b) 

where I is the identity matrix and we denote the Laplace 
transform of g(u) as g(u) = J a e~ UT g(r) dr. However, 

even though Eq. I|9a|l determines T(u) from K(u), the 
inverse Laplace transform needed to determine T(r) can 
be difficult, if not impossible, to compute. Thus, even 
though Eqs. (0 and JSJ arc formally equivalent, only 
the CTRW formulation provides a practical way to use 
mesoscopic numerical simulations to compute <p(r). 

Faradjian and Elber j^] have recently demonstrated the 
feasibility of integrating Eqs. JHJ with K(t) determined 
by molecular dynamics to compute transitions along a 
single reaction-coordinate. However, computing ip(r) by 
this approach is computationally wasteful since it is de- 
termined to high temporal resolution even though the 
experimentally relevant information is usually contained 
in only a few of its low-order moments. The unnecessary 
price paid is that the complete functional form of K(t) 
must be determined by many (expensive) numerical sim- 
ulations. 

In Sec. [H]we present a fundamental new relationship 
between the FPTD moments and those of K(t) or T(t). 
We show in Sec. IIIII that this relationship can be un- 
derstood as an extension of the steady-state flux-over- 
population method 0, of computing rate constants 
to the case of steady-state relaxation. In Sec. IIVI we 
show that combining this relationship with statistical er- 
ror analysis yields a more accurate and efficient computa- 
tional algorithm. In Sec.^we demonstrate how ip(r) can 
be modelled using a few of its moments and a Gamma 
distribution or a rational function approximation to its 
Laplace transform. 



II. THE FPTD MOMENT GENERATING 
FUNCTION 

The Laplace transform of <p(r) gives the FPTD mo- 
ment generating function: 



k=0 



(T k ))a k 



k\ 



tp(-a) 



(10) 



We assume, as is true in most cases of interest, that 
r s ' s (r) (V s, s') and ip(r) decays faster than e _QmaxT as 
r — > oo for some positive a max (corresponding to the 
slowest process in the system): 

lim r s / s (r)e Q — T = lim ip( T )e a ™* T = . (11) 
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Thus, (p{— a) is analytic in a neighborhood about and 
can be differentiated to yield all moments. (This assump- 
tion is not essential, but simplifies the discussion. If it 
is not true, the analysis will still be valid for the finite 
moments.) 

Taking the Laplace transform of Eq. JQ) gives 



ip(— a) = —aif ■ P{—a) . 



(12) 



Eq. (J2J implies that e/ • P(— a) is not analytic at a = 0, 
so expanding (p{— a) using this form is delicate |2^|. To 
avoid this inconvenience, we use Eq. to rewrite Eq. 
Q in a form that does not explicitly involve Pf(t): 



6(r- 0+) - 1 • n • dP{t)/dt\ t=T 

S(t _o+)-1- dP(t)/dt\ t=r , (13) 



where 



U = I -e 



f 



f 



is the projection operator into the dynamic subspace of 
non-absorbing states s^ / and we use the notation 

A = n-A 
m = n-Af-n 

to denote projected vectors A and matrices M. Eq. i(T3|) 
relates the FPTD to the loss of probability from the dy- 
namic states. Its Laplace transform is 



ifi(-a) = 1 + al ■ P{~a) 



(14) 



This form is advantageous because Eqs. (Q, ©, Q, and 
(|TT]l imply that 



lim P(t) e Q " 

t— >OC 



o, 



(15) 



so P(— a) is analytic at a = 0. 

To complete the solution, we express P in terms of F 
by projecting Eq. © [using U-T-P = L P, which follows 
from Eq. JSJ] and taking its Laplace transform to get 



uP(u) = Si - T(u) ■ P(u) 



The solution is 



P{u) = [uU + T(u) 



(16) 



where the use of the matrix pseudo-inverse (i.e., the in- 
verse within the dynamic subspace) is implied here and 
below. Combining this with Eq. (|14f) gives 



(p{— a) = 1 — al ■ [all — T(— a) 



(17) 



The right-hand-side is analytic at a = because F(0) is 
invertible within the dynamic subspace. 

Using Eq. 19a|l. we reexpress this in terms of K: 



(18) 



<p(-a) = -1 ■ K{-a) ■ [n + Kp(-a) 



Since K(—a) is analytic at a = even without projec- 
tion, Eq. I|18H can equivalently be written in unprojected 
form as [24| 

(p{-a) =€/•[/+ Kpi-a)]^ 1 ■ e t . 

Equations l|17() and l|18l) provide the fundamental rela- 
tionship between the FPT moment generating function 
and F and K. 



III. THE MOMENT GENERATING FUNCTION 
AND STEADY-STATE RELAXATION 

To elucidate the physical significance of Eq. (|17f) . we 
compare the computation of the MFPT using the gen- 
eralized master equation with the computation of the 
transition rate (MFPT _ ^_) using the steady-state flux- 
over-population method 0, H2 • The latter computes the 
rate as the magnitude of the flux of systems divided by 
the total dynamic population in a steady-state situation. 

We begin to relate the generalized master equation to 
the steady-state by noting that the solution of Eq. 10, 
P(t), gives the response of a linear system to an impulse 
and so is the Green's function for the general solution: 
If systems are injected continuously at a non-negative 
rate r(t) beginning at t — 0, the resultant population 
distribution vector P[r; t] will satisfy 



dP[r:t] 
Jt 



= 0{t)r{t)ii 



r(r) • P[r;t - r]dr (19) 



with boundary condition f[r, — oo] = (8 is the Heavi- 
side step- function). This has the solution 



P[r-t] 



P(t - t')r(t')dt' 



(20) 



Unlike the Green's function P(t), which satisfies 1 • 
P(t) = 1 (t > 0), the general solution P[r; t] isunnormal- 
ized; the total population is 1 • P[r;t] = 8(t) J Q r(t')dt', 
which can increase without bound as t — > oo because of 
the accumulation of systems in /. To avoid this compli- 
cation, we follow the approach used above and focus on 
the projected dynamic population vector P[r;t], which 
is bounded. 

The steady-state case corresponds to the asymptotic 
(t 3> a~ax) regime with r(t) = j, where j is a positive 
constant. More generally, we consider steady-state re- 
laxation: the asymptotic regime with r(t) — j exp(— at) 
(a > 0). Eqs. I|15f) and (|20() imply the asymptotic form 



P\je- at ;t] ~ je 



-at ■ 



.a 



(a > 0,t > aj ax ) , 



where .Poo [a] is a vector constant. Substituting this into 
the projected form of Eq. (|19fl . multiplying by exp(a£)/j 
and taking the limit t — > oo gives 



-aPoola] = e l 



e aT t{T)dr ■ P x [a] 
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with solution 

P 00 [a\=-[all-f(-a)]- 1 -e i , 

where the pseudo-inverse is again implied. Comparing 
this with Eq. (|16fl implies that 

P(-a)=P 00 [a], 

so Eq. 1)14(1 implies that 

<p(-a) = 1 + al-Poota] . (21) 

We see that the Laplace transform of Pit), and hence 
the FPTD generating function, is simply related to the 
steady-state relaxation dynamic population vector. The 
steady-state calculation of the transition rate is a special 
case of this more general relationship: Eqs. 1(10(1 and 121(1 
imply that 

«t»=1-Poo[0]. (22) 

Since jPoo[0] is the steady-state solution for constant flux 
j and its inner product with 1 is the sum over the pop- 
ulation in all the dynamic states, Eq. ((22(1 states that 
the MFPT is the total dynamic population over the flux. 
This is equivalent to the statement 0, that the tran- 
sition rate is the flux-over- (dynamic) population. 

IV. EFFICIENT CALCULATION OF THE FPTD 
MOMENTS 

Eqs. 1(11) (I and 117(1 imply that (ip) = 1, in agreement 
with Eq. (J5J. Expanding Eq. ((17(1 to first order gives 

«r)) = l.P(0)- 1 .6i = l-<f)- 1 .e i . (23) 

The same result would be obtained if we ignored all mem- 
ory effects and approximated T(r) sw 8(t)(T), which is 
equivalent to replacing the generalized master equation 
with a regular master equation having T = (V). Differ- 
ences between the moments of these two equations only 
appear in higher order. 

When only K(t) is known, we can use Eq. I(9a(l or 
expand Eq. 1(18(1 to reexpress Eq. ((23(1 in terms of K: 

{{r)) = l-{rK D )-{K)- l -e l . (24) 

To use this relationship to compute ((r)) from numerical 
simulation data, the time-averages on its right-hand-side 
can be approximated by 

n 3 

(tK d ) ss » n^gr/, (25a) 

z=l 

(K) B > S w n s > s /n s (25b) 

where n s is the number of simulations that were ini- 
tiated in state s, n s i s is the number of those simula- 
tions that made their first transition to state s' and 



{if : i — 1, . . . ,n s } is the set of first transition times 
for the simulations initiated in s. This result is much 
easier to compute than numerically solving the CTRW 
Eqs. (JHl and then integrating Eq. (|TJ to compute ((r)). 
Moreover it does not introduce quantization error, as 
occurs when numerically solving Eqs. (JSJ); its estimate 
for ((t)) equals that which would be obtained using the 
CTRW equations in the limit where the numerical quan- 
tization size h — > 0. To prove this, note that we have 
already proved that Eq. and the CTRW calcula- 

tion are equivalent when the exact K(i) and moments 
are used. These exact values would be obtained in the 
limit of an infinite amount of simulation data. The re- 
sult obtained with a finite amount of simulation data 
can be viewed as an approximation to the exact result. 
Alternatively, it can be viewed as the exact result for 
the problem in which K s ' s (r) is proportional to a sum 
of 5-functions, each corresponding to one of the wait- 
ing times in the set {if s ■ i = 1, • • • , n s i s } of numeri- 
cally computed local first passage times from s to s' i.e., 
K s , s {t) = -nj 1 Y%={ S[t -rf s ]. The numerical K(t) 
computed in the limit h — > and the moments of K 
computed using Eqs. ((25(1 are both exact for this modi- 
fied problem, and thus must yield the same result. 

A. Improving accuracy and efficiency by sampling 
adjustment 

The accuracy of the steady-state computation of ((r)) 
will depend on the quality of the statistical estimation of 
(tKd) and (K) provided by Eqs. The simplest way 
to use these equations would be to follow the procedure 
used to estimate K{t) in the CTRW approach and to 
initiate the same number of simulations in each state; 
i.e., n s — n to t/{N — 1), where n to t is the total number of 
simulations to be performed. (The denominator is N — 1 
because no simulations are initiated in the final state.) 
However, this procedure is not optimal because it does 
not account for differences in the sensitivity of the result 
to errors in different states. For example, the inverse 
matrix (K) -1 appearing in Eq. (|2"4l can be particularly 
sensitive to errors in small matrix elements corresponding 
to bottlenecks in the probabilistic flow where there tend 
to be fewer transitions in the "forward" direction. Since 
the expected root-mean-square (rms) statistical errors of 
the matrix elements {K) s i s are inversely proportional to 
n s , overall accuracy will be improved if n s is increased 
for the bottleneck states while being decreased for other 
states to keep n to t constant. 

We can use Eq. 1(24(1 to analyze the dependence of a 2 , 
the variance of ((r)), on the n s and thereby to quantita- 
tively optimize effort allocation. To simplify notation we 
define 

p s s , = -(Ks's), 

4> s = {tK d ) ss . 

The p s s , are multinomial probabilities governing first tran- 
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sitions out of state s, which by Eqs. (J5J) and J7J) satisfy 
J2 s 'yisPs' ~ 1- Accounting for the reduction in the stan- 
dard error of the mean resulting from repeated sampling, 
making the approximation that the statistical errors in 
(j) s are independent of those in the p s s , |25| , and using the 
propagation of errors formula, we estimate 



- 2 = E^ 



d((r}) 



E 



a«7 



-a* 



d((i 



where a 2 . s is the variance of the {r/ } and 



dp s s , 
(26) 



-Pi 



(s" = s' 
(s" ± s> 



is the approximate multinomial variance tensor for state 
s |26|. Since the cost of a simulation is proportional to 
its duration, the expected cost of n s simulations initiated 
at state s will be n s (j) s . Minimizing a 2 with respect to 
the n s while maintaining a constant total cost implies the 
optimality conditions 



'1/2 



\ 



9{{r)) 



E 



d{{n 



9«7 



dp s s „ 



s;s " s ' dp s s , 

where c is a constant chosen so that S s ^/ n s4>s = cost. 
Eq. 1)27(1 determines the n s as explicit functions of the <f> s , 
<J 2 . S , and p s s , . To estimate these parameters, we can first 
perform a pilot run with a small number of simulations 
for each state. More simulations can then be added to 
the pilot simulations so that the combined set satisfies 
Eqs. (|77jl ^23. Eq. can be used to estimate the error 
of the final result computed with the combined set of 
simulations to determine if the accuracy goal has been 
met. 

Empirically, we have found that efficiency can be fur- 
ther improved to a small extent by replacing the max- 
imum likelihood estimator of the p s s , used in Eq. I|25b|) 
by a Bayes-Laplace estimator (Appendix A). This esti- 
mator was used in the example discussed below but only 
gave noticeable improvement for the low-accuracy (e.g., 
25-50%) results jig. 



B. Example 

We compared the efficiency of the sampling-adjusted 
steady-state procedure with that of the standard CTRW 
procedure using the two-dimensional entropic barrier 
model studied by Faradjian and Elber [Sj- They com- 
puted the FPTD for transitions under Brownian dy- 
namics with potential energy function U(x, y) — x e + 
y 6 exp(— 10Cte 2 )[l — exp(— 100z/ 2 )] from an initial state 
with x = —1 to a final state with x = 0.714 at kT— 
0.5 and friction coefficient 7 = 0.1 The "exact" value of 
((t)) was computed using the CTRW method with five 



linearly-ordered intermediate states and n s = 5,000 nu- 
merical simulations initiated at each state (i.e., a total 
of Utot — 30, 000 simulations were used with cost = 
maxcost = 5, 000 Y^s^f ^ s )- To assess the accuracy of 
the method as cost was decreased, we used their simu- 
lation data to determine the geometric rms error |30| of 
the CTRW estimates of ((r)) when fewer simulations were 
used corresponding to maxcost/cost = 2, 4, 8, 16, 32, 64, 
128, and 256. For each value of cost CTRW estimates of 
((t)) were computed for 4,000 random data subsets, and 
their geometric rms error was computed relative to the 
"exact" value [29j. 

To assess the performance of the sample-adjusted 
steady-state procedure for a specified cost, we first per- 
formed a pilot run (with n s the same for all states) 
costing 1/4 cost, used the estimated values of <fi s , cr 2 ^, 
and p s s , and Eqs. I|27|l to optimize the distribution across 
states of additional simulations costing 3/4 cost, and 
evaluated ((r)) using Eqs. ||2U), l|25a|) . and the Bayes- 
Laplace proportions estimator (Appendix A) with the 
combined set of simulations. This procedure was re- 
peated 4,000 times to estimate the geometric rms error. 
Additional tests showed that the results were not highly 
sensitive to the size of the pilot run. 

The geometric rms errors for both methods as a func- 
tion of cost are plotted in Fig. 1 and show that the sam- 
pling adjustment increased efficiency slightly more than 
two-fold. For example, cost = maxcost/8 was needed 
to achieve ~ 12% accuracy using the standard CTRW 
method, while only maxcost/16 was needed for ~ 11% 
accuracy with the adjusted steady-state method. Exami- 
nation of the optimized n s showed that this gain occurred 
because a ~ 4-fold increase in the sampling frequency at 
a bottleneck caused a ~ -\/4-fold improvement in the as- 
sociated dominating error. 

The extent of sampling adjustment in this problem was 
limited because there were only five mesoscopic states 
among which effort could be reallocated. Larger adjust- 
ments, and larger gains in efficiency, may be possible in 
larger problems if the increase in the number of meso- 
scopic states exceeds the relative increase in the number 
of bottlenecks. Such gains could be particularly impor- 
tant for very costly problems (e.g., those arising when 
studying protein conformational transitions). 

Expressions analogous to Eq. (|24|l for the higher FPT 
moments can be obtained by analytically expanding Eq. 
(fTSfl in terms of the (r k K) .31] ■ Although the optimal 
sampling conditions for simultaneously computing multi- 
ple moments differ from Eqs. I|27|l . we expect that benefit 
will still be achieved even if sampling is adjusted using 
these equations. Of course, even better results will be 
obtained if the optimization analysis is extended to the 
multiple moment case. 
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V. EFFICIENT MODELLING OF THE FPTD 

In some cases we will need to compute not just the mo- 
ments, but also ip(r) to a low temporal resolution com- 
mensurate with experimental results. We can extend the 
efficiency improvement obtained in the moment compu- 
tations to this case by modelling the FPTD using a few 
of its low-order moments and an appropriate functional 
expansion. 



A. Modelling using the Gamma distribution 



A Gamma distribution of the form 

f(P,r,T 



/ 3(/3r)''- 1 e-' 3T 



r(7) 

(where (3 > 0, 7 > 0; here T denotes the Euler Gamma 
function, not the transition matrix) provides a simple 
model. It decays exponentially as t — > 00, thereby match- 
ing the expected asymptotic behavior of <p(t), and it is 
simple to choose (3 and 7 so as to match the first two 
moments of the FPTD 



(r k f) = 



k = 1,2 



(T /) - T 



13 = 7/«t» 



(28a) 

(28b) 
(28c) 



[The denominator of the expression for 7 is equal to (((r— 
((t))) 2 )), and so is guaranteed to be positive.] Additional 
moments could be included by modelling ip(r) as a sum 
of Gamma distributions, but problems with non-unique 
parameter fitting can arise. 



B. Modelling using a rational function 
approximation to (p 

In some cases better results can be obtained by ap- 
proximating (p{u) as a rational function 



<p(u) « i? m> „ 



(u) = 



1 



Piu - 



■ p m u" 



1 + qm + 



■ q n u" 



(n > m) 



(29) 

Here we have fixed the zeroth-order terms in the numer- 
ator and denominator so that R m ,n(0) = <f(0) — 1, as 
required by Eq. (J2J. We require n — m > 1 so that 
lim u _ >oo R m . n (u) vanishes at least as fast as s -2 , implying 
that its inverse Laplace transform will vanish at the ori- 
gin, corresponding to f(0) = 0. Since the only singulari- 
ties of R m ,n(u) are poles, the inverse Laplace transform 
is easy to compute. Moreover, if the only poles are on the 
negative real axis (not guaranteed) , the inverse transform 
will be the sum of decaying exponentials, thereby provid- 
ing a natural model for tp(r). We use this property as a 



validity check and do not accept (potentially overfitted) 
approximations that have poles off the negative real axis. 

The pi and qi are fixed by requiring that <p(u) = 
Rm.n {u) at m+n non-zero values of [k = 1, . . . , m+n) . 
To choose the {ufe} appropriately we note that the most 
important structure of ip(r) occurs at scale t ~ ((■?"))■ 
Therefore, the important structure of (p{u) will occur at 
scale u ~ ((t))~ . Thus we choose Uk = k/((r)) and re- 
quire 

-fcr/«T»u 



<p(k/((T))) = {(e 



flm >B (*/H) 

(k = 1, . . . , m 



(30) 



■n). 



The statistical error of the exponential moments grows 
as k increases because the exponential will down-weight 
a larger fraction of the data points. This limits the accu- 
racy of the "high-frequency" components of the moment- 
modelled <p{r) to be the same as that of the directly in- 
tegrated ip(r) [33|]. 

The first moment of (f obtained using Eq. (|29|l will be 
close to, but will not exactly match the MFPT. An exact 
match can be obtained by replacing Eq. (|29|l with the 
constrained rational function 



(p{u) » Rthu) 



1 + piu 



■PmU" 



1 + (((r)) + pi)u + q 2 u 2 + ... q n u n 

(n > m) 



(31) 



This satisfies 



da 



(32) 



so the first moment of its inverse Laplace transform will 
exactly equal the MFPT. The constraint on the deriva- 



in Eq. 

At)) 



At)) 

(JSOJl, so when Eq. 



tive of replaces the use of the k = 1 constraint 



Rm,n(uk) = (p{Uk) for fc 



l|31|l is used we only match 
= 2, . . . ,m + n 34]. In most 



cases the estimates obtained using Eqs. or Eqs. lJ3T|) 
will be similar. The two lowest-order approximations of 
this type are those involving the MFPT and either one 
or two exponential moments corresponding to approxi- 



mating ip (a) as i?Q2^(a) or R^'^a). 

We illustrate the method using the two-dimensional 
model discussed above. In Fig. 2 we compare the FPTD 
computed by Faradjian and Elber using the CTRW 
method with the approximated FPTD's computed us- 
ing the MFPT and either zero, one or two exponential 
moments. The single-exponential fit obtained using the 
MFPT alone [i.e., exp(— f/((r)))] misses much important 
detail, but a fairly good representation is obtained by ad- 
ditionally matching just one exponential moment using 
ip(a) « R { $ } {a). The fit obtained with the MFPT plus 
two exponential moments (m = 0, n = 3) is practically 
indistinguishable from the exact ip(r). The next higher 
order approximations have imaginary poles. This pro- 
vides an indication of overfitting and (correctly) suggests 
that the approximation should not be extended further. 



At)), 
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While this procedure has worked on a few tested cases, 
as with all parameterized modelling approaches, success 
depends on a reasonable match between the form of the 
parameterized approximation and the true distribution. 
This can not be guaranteed but is a reasonable assump- 
tion since most FPTDs are expected to have distributions 
qualitatively like that shown in Fig. 2. 

VI. SUMMARY 

Mesoscopic coarse graining and the CTRW equations 
can be used to compute the macroscopic FPTD, ip{r), of 
a complex stochastic system from short-term, and hence 
affordable, microscopic numerical simulations of its dy- 
namics. In many cases interest will focus on the MFPT 
and possibly a few additional low-order FPTD moments. 
Instead of integrating the CTRW equations over time 
to compute <p(t), a procedure that requires the full func- 
tional form of K(t) to be estimated, and then integrating 
again to compute the moments, we have shown that the 
FPTD moments can be computed simply and directly 
from the moments of K . This method is simpler and 
eliminates the quantization error inherent in the numeri- 
cal solution of the CTRW equations in the time-domain. 
It can physically be viewed as an adaptation and exten- 
sion of the steady-state flux-over-population method of 
computing transition rates, so we call it steady-state re- 
laxation. 

The steady-state expressions for the FPTD moments 
are simple enough for straightforward statistical error 
analysis, which permits the accuracy of the computed 
moments for a given amount of simulation data to be 
estimated. This analysis can also be used to optimize 
the allocation of computational effort over the different 
mesoscopic states and to thereby reduce the total cost 
of the numerical simulations required for fixed accuracy. 
This is important since computability of the FPTD in 
large problems will often be limited by this cost. Such 
optimization improved efficiency over two-fold in a test 



problem with five mesoscopic states, and greater im- 
provements are possible in problems with more states. 
This improvement can be extended to the FPTD itself 
by modelling it using either a Gamma distribution or a 
rational-function approximation to its Laplace transform. 
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APPENDIX A: BAYES-LAPLACE ESTIMATOR 

In some cases (e.g., when numerical simulations are 
particularly costly), the goal may be just to estimate ((t)) 
to rough accuracy (e.g, 25-50%) using the smallest pos- 
sible number of simulations. In such cases the n s may 
be small and there may be large fractional errors in the 
(K) s ' s whose effects are amplified by the matrix inver- 
sion. Because the inversion is nonlinear, the maximum 
likelihood estimator used in Eq. I|25b|) may not be opti- 
mal and it is worth considering other possibilities. One 
alternative is the Bayes-Laplace estimator [35j j 



where v s is the number of states to which s can make 
transitions. This is the mean Bayesian estimate of (K) s i s 
using a non-informative prior distribution (i.e., making 
the a priori assumption that a system in state s is equally 
likely to make a transition to any of the connected states 
s'). This estimator has a bias away from very small 
(K)s's, suggesting that it may reduce the error of the 
inverted matrix. This surmise was empirically found to 
be true in the example of Sec. II V 51 but the improvement 
was only noticeable when the error was > 25% [28) . 
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However, this requires estimation of the (f> B i s and the 
variances of the {rf s }, &^.. s i B , which will require larger 
amounts of simulation data than estimation of <f> s and 
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timate the additional parameters, but it could be less 
effective than Eq. I|27^ for optimizing inexpensive low- 
accuracy results if not enough simulations were available 
to accurately estimate them. This is probably splitting 
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FIGURES 
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FIG. 1: Dependence of the geometric rms error of the calcu- 
lated MFPT on the total cost of the simulation for the CTRW 
and sampling-adjusted steady-state methods applied to the 
two-dimensional entropic barrier model studied by Faradjian 
and Elber The error for the CTRW method can not be 
accurately estimated for n s = 39 and n s = 19 because too 
many trials have sparse patterns of transitions that fail tjo 
connect the initial and final states* 
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FIG. 2: Comparison of the approximated and exact FPTDs 
for the two-dimensional entropic barrier model studied by 
Faradjian and Elber |<J. The exact FPTD is displayed (solid 
line) along with the approximated FPTDs obtained by ap- 
proximating the Laplace transform by a rational function 
matched to the MFPT alone (fine dotted line) or to the MFPT 
plus one (dashed line) or two (dotted line, indistinguishable 
from solid line) exponential moments. 



